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We study the saturation near threshold of the axisymmetric magnetorotational instability (MRI) 
of a viscous, resistive, incompressible fluid in a thin-gap Taylor-Couette configuration. A vertical 
magnetic field, Keplerian shear and no-slip, conducting radial boundary conditions are adopted. The 
' weakly non-linear theory leads to a real Ginzburg-Landau equation for the disturbance amplitude, 

, like in our previous idealized analysis. For small magnetic Prandtl number {Vm 1), the saturation 

' amplitude scales as Vm'^ while the magnitude of angular momentum transport scales as Vm^ ■ The 

difference with the previous scalings (oc "Pn/^ and Vm respectively) is attributed to the emergence 
of radial boundary layers. Away from those, steady-state non-linear saturation is achieved through 
a modest reduction in the destabilizing shear. These results will be useful to understand MRI 
laboratory experiments and associated numerical simulations. 

PACS numbers: Valid PACS appear here 

(N 
> 

(N 
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, The magneto-rotational instability (MRI) is a linear instability known to occur in rotating hydromagnetic shear 

■ flows when the angular velocity decreases with distance from the rotation axis, i.e. (3fl(il^) < 0. Although it had 
been known for almost half a century P, [3, Q , the MRI acquired a renewed interest only after the influential work of 
Balbus & Hawley fJ], who have shown, by means of linear stability analysis and numerical simulations, its viability 
in conditions locally approximating astrophysical accretion disks. Subsequent investigations of this kind (see the 

2 , reviews by Balbus & Hawley @ and references therein) have quite convincingly demonstrated that this instability 

■ can drive magnetohydrodynamical (MHD) turbulence in a variety of conditions, appropriate to accretion disks and 
, more general settings as well. Within the framework of a magnetic Taylor-Couette configuration, which is relevant 
' for the present work, the parameter dependencies (magnetic Reynolds and Prandtl numbers) of the marginal (linear) 

. 1^ ' instability threshold has been previously considered It was found, among other things, that the critical magnetic 

Reynolds number does not scale with the magnetic Prandtl number, for small values of the latter. This result carries 
^ ' over into the weakly nonlinear theory presented here. 

Accretion disks are important and ubiquitous astrophysical objects and are thought to power as diverse systems 
as young stellar objects, close binary systems and active galactic nuclei. Accretion disks are flattened, high specific 
angular momentum (with essentially a Keplerian distribution) masses of gas, through which matter accretes onto a 
central object. An efficient dissipation and transport of angular momentum mechanism is needed in order to allow 
accretion and reconcile theoretical models with observations. Since the typical hydrodynamical Reynolds numbers (TZ) 
in these astrophysical flows are enormous, it has been recognized at the outset, when accretion disks were theoretically 
proposed [9, 10], that some anomalous, enhanced (conceivably turbulent) dissipation and transport must be invoked. 
Keplerian rotating flows are (according to the Rayleigh and other criteria) linearly stable and thus astrophysical 
disk turbulence can not originate from a linear instability of the kind known (and well studied) in Taylor-Couette 
hydrodynamical flows. 

The physics of the non-linear development of the MRI, its saturation and the nature of the resulting angular 
momentum transport are quite complicated. Almost all of our present knowledge on this subject comes from numerical 
simulations, carried out by several groups (see, e.g., 0, [E H, and references therein). These finite-difference 
simulations, even though intended for the study of the MRI in its astrophysical setting, were actually local, i.e. done 
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for a small portion of an accretion disk, in what is known as the shearing box or sheet (hereafter, SB) formulation [ll|, 
Q (see the Appendix of [l^ for a formal account on this approximation). Although a lot has been learned from these 
simulations, the intricate processes at work are not yet fully understood and some basic physical questions remain 
open (see, e.g., [T5l|). As a result, there has recently been a growing interest to observe the instability in the laboratory, 
where various physical aspects can be unraveled in a controlled way. A number of groups have indeed embarked on 
such experimental projects, in several setups, often accompanying them by appropriate numerical calculations (e.g., 
[H, m in, m, [13, El , and references therein) . 

In comparison to the large extent of numerical and experimental work on the MRI's nonlinear development, there 
have only been very few reports on analytical and semi-analytical studies on this subject. This fact seems surprising, 
because a very large body of work, utilizing various asymptotic approaches, has been done for other important fluid 
instabilities (for reviews see, e.g., [H, [H, [H, [H, H^). We are aware of only two asymptotic studies of this kind in 
the MRI context: 

• In the first one [131 , Knobloch & Julien investigated the saturation of the MRI in the strongly nonlinear (far 
from instability threshold) regime. They utilized the so-called channel modes (radially independent axisymmetric 
linear modes, which also happen to be exact solutions of the nonlinear problem in the SB formulation [5j.[28j). 
They performed an asymptotic calculation, in which the evolution of channel modes is followed into the nonlinear 
regime by gently tuning the system out of the developed short-wavelength channel mode configuration (and 
under a specific regime of system's parameters). This work shows that nonlinearities saturate the system in 
such a way that the momentum transport scales as {TZTZm)~^ , where TZ and TZm are the hydrodynamic and 
magnetic Reynolds numbers, respectively (see their Eq. 4.22). The results further indicate that, by modifying 
the underlying shear (the "source" of instability), the system saturates while approaching solid body rotation. 

• In the second study [2^ . hereafter UMR06, we have employed a more traditional approach - weakly nonlinear 
asymptotics close to the instability threshold. The problem we considered differed from previous studies in that 
we considered the dynamics to be restricted to a narrow (in its radial extent) channel. Our original intent was 
to understand the MRI under a more controlled setting - one in which the channel modes are filtered out by the 
imposition of no normal-flow conditions at the inner and outer boundaries of the channel. Under these conditions, 
arguably more appropriate to capture the physics of experimental setups, the MRI unstable mode transits into 
instability in a way analogous to that of Rayleigh-Benard convection. An idealization, involving a hybrid 
free-slip/no-slip and conducting/insulating boundary conditions, atop the no-normal flow conditions mentioned 
above, allows for transparent analytical evaluations of the derived necessary quantities (similar idealizations have 
sometimes been used in other studies [20| ) of the problem. The similarity of this formulation to other extensively 
studied hydrodynamical instability problems led us to the application of weakly nonlinear asymptotic techniques 
to examine the system's transition into the nonlinear realm, as well as to comparison of the results to specially- 
designed numerical simulations. We found that, as the system is gently tuned into instability (through a suitably 
defined non-dimensional parameter e) , a saturated pattern-state emerges with an amplitude of the most unstable 
mode evolving according to the real Ginzburg-Landau equation (GLE), 

d^A^\A + DdlA-a\A\^A^Q, (1) 

where T and Z are suitably "stretched" time and vertical coordinates and the coefficients of the equation are 
all real and computable from the parameters of the physical problem. In particular, the coefficient a was found 
to scale as V^, where Vm is the magnetic Prandtl number, defined by the ratio TZ„i/TZ. It means that the 
amplitude achieved by the system in the saturated state scales as e^JVm and correspondingly, the overall angular 
momentum transport as e^Vm- For TZm fixed this transport would scale like l/TZ and this formulation is useful 
when the resistivity of the medium is set by its physical state (i.e. degree of ionization) and one wishes to 
estimate the effect of decreasing effective viscosity (resulting, e.g. from the inaccuracy of the numerical scheme 
in a simulation). These analytical scalings were found in the limit Vm «C 1, while for larger values oiVm similar 
trends may be expected but the coefficients have to be evaluated numerically. We have conjectured that for 
self-consistent boundary conditions the above general qualitative behavior should hold as well, with perhaps 
some change in the relevant power of Vm in the scalings. Our asymptotic analysis was accompanied by fully 
numerical spectral calculations of the original SB equations with similarly idealized boundary conditions. The 
analytical and numerical scalings were found to agree quite well. 

In this paper we present a study of the MRI as developing in a model representing the thin-gap limit of a magnetic 
Taylor- Couette (hereafter mTC) configuration, in which an incompressible axisymmetric rotating flow is subject to an 
external vertical magnetic field. This will permit a quantitative examination of the effect of the boundary conditions 
on the results reported in UMR06 and confirm the conjecture on the general qualitative behavior. 
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The fundamental equations of motion are the same as those assumed in previous studies of the MRI (e.g. [4|) 
save for the inclusion of non-ideal effects, namely resistivity and viscosity. Solutions to these equations are sought, 
subject to realistic boundary conditions at the system walls, namely that of no-flow and conducting conditions. For 
the vertical boundary conditions we assume periodicity for the sake of simplicity and transparency. After presenting, 
in Section ini the relevant approximations, definitions and equations, we perform, in Section [llll a linear eigenmode 
analysis. We identify the most unstable mode as a function of the non-dimensional parameters of the system - of 
which there are five: the Cowling number C, the magnetic Prandtl number Vm, the magnetic Reynolds number TZrm 
and shear index q (see below). We demonstrate next that this system has a transition into instability which is similar 
in some important aspects to that in Rayleigh-Benard convection [l^, [2^ [2^ . We also identify the presence of a 
neutral, spatially constant mode representing the hand of a constant azimuthal field. 

In Section IIVI we perform a weakly nonlinear asymptotic analysis by tuning the system away from the conditions 
of marginality. In this case this is done by ratcheting the background magnetic field downward from the marginal 
state with the magnitude of the departure from that state measured by the small parameter e^. The full calculation, 
detailed in Appendices B-D, reveals that the envelope (of the marginally unstable modes) evolution is governed by two 
uncoupled partial differential equations: one represents the leading MRI mode and evolves according to the real GLE 
and the other equation, representing the evolution of the uniform azimuthal field, is a standard diffusion equation. 

2/3 

The saturated amplitude of the leading MRI mode is demonstrated, in the Vm <C 1 limit, to scale as eVm and is 
shown to be affected by the boundary layers appearing at the system walls. The main physical factor contributing 
to saturation is identified as coming from the second order (in e) correction to the azimuthal velocity perturbation in 
the limit Vm ^ 1- This, in turn, affects the shear profile so as to stabilize the new steady configuration. We also find 
that the average total angular momentum transport implied under these conditions scales as e^Vm^ for Vm <C 1, or 
as t^TZ~^/^ for TZm fixed (and of O (1)). These results are in accord with our conjecture and expectations given in 
UMR06. 

In the last Section we discuss the implications of our work and how it should be perceived as a part of the ongoing 
research efforts on various aspects of the MRI. Wc also provide some heuristic arguments to help understand the 
results. Finally, we end with a short outline of possible directions for future work of this kind. 



II. ASSUMPTIONS, DEFINITIONS AND EQUATIONS 

The hydromagnetic equations in cylindrical coordinates fS] are applied to the neighborhood of a representative 
radial point (rg) in the system, using the above mentioned shearing box (SB) approximation. The SB is applied here 
to the thin-gap limit of a Taylor-Couette setup with an imposed background vertical magnetic field. We begin by 
considering a steady base flow with only a constant vertical magnetic field, B — BqZ, and a velocity of the form 
V = U{x)y. In this base state the velocity has a linear shear profile U{x) = —qftox, representing an azimuthal fiow 
about a point tq, that rotates with a rate flo, defined from the differential rotation law fi(r) oc flo^r /ro)~'^ ■ The total 
pressure in the base state (divided by the constant density). 



Po \ Stt 

is a constant and thus its gradient is zero. 

This base flow is disturbed by 3-D perturbations on the magnetic field b = {bx,by,bz), as well as on the velocity - 
u = (ux,Uy,Uz), and on the total pressure - ■uj. We consider only axisymmetric disturbances, i.e. perturbations with 
structure only in the x and z directions. This results, after non-dimensionalization, in the following set of non-linear 
equations: 

— - 2noz X u - qflou^y - Cb • Vb - CB^d^h = -Vw + — V^u, (2) 
dt TZ 

^~h-yu + qnob,y-Bod,u - :5^V2b, (3) 
dt Km 

together with an incompressibility condition and the solenoidal magnetic field constraint 

V • u = d^ux + d,u, = 0, V • b EE d^bx + d,b, = 0. (4) 

The Cartesian coordinates x, y, z represent here the radial (shear-wise), azimuthal (stream- wise) and vertical directions 
respectively and since axisymmetry is assumed V = yidx + zdz and the Laplacian is = 9^ -I- 9^ . Lengths have been 
non-dimensionalized by L (the shearing-box size), time t by the local rotation rate fJo (tildes denote here dimensional 
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quantities). Because the dimensional rotation rate of the box (about the central object) is CIq = ttoz, the non- 
dimensional quantity f2o is formally equivalent to 1, but we keep it to flag the Coriolis terms. Velocities have been 
scaled by ripL and the magnetic field by the value of the background vertical field Bq. Thus the non-dimensional 
constant background field Bq = 1, but again, we leave it in the equation set for later convenience (see below). The 
hydrodynamic pressure is scaled by poL^ilg and the magnetic one by _BQ/(87r). The non-dimensional perturbation w 
of the total pressure divided by the density (which is equal to 1 in non-dimensional units) , which survives the spatial 
derivatives, is thus given by 



w = p - 



Ci|bp, (5) 



where p is the hydrodynamic pressure perturbation. 
The non-dimensional parameter 



^Xl (6) 

is the Cowling number, measuring the relative importance of the magnetic pressure to the hydrodynamical one. It 
is equal to the inverse square of the typical Alfven number {Va is the typical Alfven speed). The Cowling number 
appears in the non-linear equations, together with the two Reynolds numbers 

n = ^^, 7e™ = ^— , 7 

1/ r] 

where v and fj are, respectively, the microscopic viscosity and magnetic resistivity of the fiuid. We shall also see 
that the magnetic Prandtl number, given as Vm = T^rn/Ti-, plays an important role in the nonlinear evolution of this 
system. 

We rewrite now the equations of motion in terms of more convenient dependent variables: 

dtV^-ii + N^ ^ n^'v^'ii + 2nod^Uy+CBod^\7^<P (8) 

dtUy + N^ ^ n\'^Uy-nQ{2- q)d^^ + CBod,^hy (9) 

dt^ + N^ = 7e„"v2$ + So9,* (10) 

dtby + N^ = n,n\^by + B^d.Uy - qO^d^, (11) 

Because the fiow is incompressible and ^/-independent, the radial and vertical velocities are expressed in terms of 
the streamfunction, that is, {ux,Uz) = {dz'^ , —dx'^)- Also, since the magnetic field is source free, we similarly 
express its vertical and radial components in terms of the flux function, <i>, that is, {bx, bz) = {dz^, —dx^)- Note that 
(1101) combines information about the radial and vertical magnetic fields in terms of the flux and streamfunctions (e.g. 
[20|). In this formulation the nonlinear advection and tension terms are 

A^^ EE J(*,V2*)-CJ($,V2$), N^= Ji-^,Uy)-CJi<P,by), 

N,= J(*,<i>), N,=j{^,by)-J{<i>,Uy), (12) 

in which the Jacobian is defined as J{f,g) = dzf dxg — dxf dzg- The underlined term in (jlip . representing the 
transport of the perturbed radial magnetic field by the background shear flow, is instrumental for the occurrence of 
the MRI in this system. 

The boundary conditions are periodic on the vertical boundaries of the domain and we require also that the flow 
be no-slip at the inner and outer boundaries. This means that u = at a; = ±1, i.e. 

Uy = 0, = 0, dx'i' = 0, at X = ±1. (13) 

Regarding the boundary conditions on the magnetic field disturbances, we posit conditions (only 2 are needed) that 
are consistent with the inner and outer walls being conducting, bx ~ and dxby = at a; = ±1, i.e. 

dz<S> = 0, dxby = 0, at X = ±1. (14) 

Note that these boundary conditions arc more physically consistent than the ones we have used in UMR06, however 
they will call for a numerical evaluation of the eigenfunctions and the coefficients for the asymptotic analysis that 
result from them. 
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FIG. 1: (Color online) Summary of linear theory. This example is for C = 0.08, Vm = 0.001, q — 3/2, and the fundamental 
mode, (a) Growth rates, Re{s), as a function of wavenumber k for three values of TZm- (b) Solid line depicts those values of TZm 
and k where Re{s) = 0. The shaded region shows unstable modes. The locations of fc = kcrit = Q and TZm = 7?.m(crit) = R,n 
are shown. 

Finally, we point out that there exists an energy theorem for the above dynamical equations. Defining the total 
energy (per unit length in the azimuthal direction) of the disturbances in the domain as £' = i J (u^ + Cb^) dxdz, 
we get, after the usual integration procedures and application of boundary conditions, 



dE 
'dt 



qflo / Tdxdz 



1 

n 



where 



T = T 



iVuzl ) dxdz - 



_C_ 

TP 



(|V6, 



R 



M, 



Tr = UxUy, Tm = -Cbxby. 



\Vb^\^)dxdz, (15) 



Tr and Tm are the Reynolds (hydrodynamic) and Maxwell stresses, capturing the velocity and magnetic field dis- 
turbance correlations, respectively. Statement (jlSp is analogous to the Reynolds-Orr relation in hydrodynamics (for 
which T — Tr and Tr — 0). The total stress T will be used in the asymptotic theory we develop here as the dominant 
expression for the evaluation of transport, occurring during the weakly nonlinear evolution of the system. The full 
RHS of including the two dissipative terms, should obviously vanish when a saturated, steady state is reached. 
We discuss this in more detail in Section V. 



III. LINEAR THEORY 



Linearization of ([51fTT|) yields the following equation 



(16) 



in which all the small perturbations are lumped in the vector Vi = (^'^ (x), (x), (x), 6^ (a;)) e + + cc, with 
k being the vertical wave-number and s the temporal eigenvalue. The spatial differential operators V and C (appro- 
priately written in the form of 4 x 4 matrices) are explicitly given in (B1-B3) of Appendix B. As long as A: 7^ the 
boundary conditions on the functions of Vi become (see eqs. []"3lfT4)) 



0, 



at x = ±1, 



(17) 



where = d/dx. 

In principle, equations (|16p can be set up and solved analytically however the resulting expressions are far too 
cumbersome to be conveniently manipulated. It is much easier to solve this set numerically, using a Chebyshev 
collocation technique. Each function is approximated using typically between 30 and 60 grid points on a Chebyshev 
numerical grid. Larger number of points are required for smaller values of the magnetic Prandtl number. 

We shall concentrate on and follow here only one mode and call it the fundamental one. This is the mode which first 
becomes unstable when the vertical magnetic field is decreased below threshold (the mode is marginal at threshold). 
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For given values of the parameters, the eigenvalue corresponding to this fundamental mode arises as one of the four 
possible solutions of the dispersion relation. It is purely real (/m(s) — 0) and thus the instability is steady, or non- 
oscillatory (in the customary nomenclature, e.g., [23|). The solution of the dispersion relation provides the functional 
dependence s — s{k,q,C,Vm,TZm)- 

In Figure ([T]-a) we display the growth rate i?e(s) as a function of k of this fundamental mode, for several values 
of TZm- The parameters C,Vm and q are fixed at the values indicated in the caption. We see that the transition 
into instability is typical of steady-cellular instabilities (similar, in principle, to Rayleigh-Benard convection). The 
marginal mode can be chosen to have a transition to instability at the maximum of the curve s{k) (i.e. s = 
simultaneously with ds/dk = 0), while all the other modes show strong temporal decay. The marginal mode can 
be identified with respect to a critical wavenumber fccrit = Q and a critical magnetic Reynolds number 7?.m(crit) = R^- 

Figure ([T]-b), which shows the neutral curve (s = 0) in the TZm — k plane, also demonstrates the way in which 
the critical values Q and Rm are determined. These critical parameters are in general functions of the remaining 
parameters of the system, i.e. Q = Q{q,Vm,C) and Rm = RmiqiVmTC). From here on out we will restrict our 
considerations to values of g = 3/2 (for consistency with UMR06) and consider the behavior of these quantities as a 
function of C and, primarily, Vm- 

The eigenfunctions for the mode in question have even symmetry with respect to x = due to both the symmetry 
in the boundary conditions and the symmetries inherent to the thin-gap limit of the mTC problem. In Figure ^ we 
display a sample of eigenfunctions of the marginal mode. To avoid later notational ambiguity, the eigenfunctions for 
these marginal modes (i.e. those with k — Q and TZm = Rm) will be labeled with a "11" subscript, that is, those 
modes will be represented by 

*i ^ <-> -Uji, $1 ^ $11, 6i ^ fell, when k = fccrit = Q, T^m = 7?.,„(crit) = Rm- 

It is argued in Appendix [Xl that in the limit T'm <C 1, the boundary layers size that appear scale as T'm^- The 
boundary layers that develop are satisfactorily represented numerically by the Chebyshev method used, e.g. with a 
grid of 50 points we can resolve at least 3-4 points of the boundary layer zones when T^m = 10^^. This dependence 
on T^m will also have some bearing on the scaling properties of the coefficients of the resulting (real GLE) envelope 
equation, presented in the next section. 

Finally, we note that there always exists an additional marginal mode of the system, separate from the above 
mentioned MRI mode. This neutral mode reflects a symmetry introduced into the system due to the conducting 
boundary conditions. Namely, a spatially constant, time-independent solution for the azimuthal magnetic field (i.e. 
by = constant) solves both the linear (and, incidentally, the nonlinear) equations and satisfies its requisite boundary 
conditions. This mode must be formally included in the subsequent nonlinear analysis. 



The weakly nonlinear analysis aims to develop a description of the system's evolution beginning very close to 
marginality, slightly into the unstable region. The control parameter in the asymptotic analysis is incorporated in 
the expression for the background magnetic field. It is here set to be i?o = 1 ^ e^, i-e. the degree of departure from 
marginality is controlled by the small parameter e (of our choosing) whose only formal restriction is that it he e <^ 1. 

Close to marginality the relevant MRI mode, discussed in the previous section, may be expressed to leading order 
in e (as can be shown by a simple scaling and balancing analysis) in the form 



where Vn = (^'n , Mh, , fen)"^ and Un = (0, 0, 0, 1)"^. The inclusion of BV^^ in this general solution is dictated by 
the presence of the neutral mode, discussed at the end of the previous Section. 

The weakly nonlinear evolution is asymptotically derived by allowing the amplitudes A and B to be (weakly) 
dependent on space and time. The aim is to develop an evolution equation for the envelopes A and B (space and time 
dependent amplitudes) as one tunes the system away from the marginal state defined above &t k = Q and TZm = Rm- 
The wisdom garnered from other problems involving cellular instabilities [12, [13, [13j [1^ guides us into an Ansatz 
such that the two envelope functions have functional dependencies upon a long time scale, T = e^t and a long vertical 
scale, Z = ez, i.e. we posit the form A = A{e^t, ez), B — B{e^t, ez). The end-result of this asymptotic procedure, fully 
detailed in Appendices B and C, are the two (decoupled) amplitude equations 



IV. WEAKLY NONLINEAR ASYMPTOTIC ANALYSIS 



eVi =e(AViie*'5"+BUii +c.c), 



d^A = XA~aA\A\^ +DdlA, 



(18) 




(19) 
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FIG. 2: (Color online) Eigenfunctions for the marginal mode at C = 0.08, P,„ = 0.001, q = 3/2. Here Q ~ 0.75, 7?„ ~ 4.9. The 
eigenfunctions are shown fitted (solid line) to the values determined numerically on the Chebyshev grid (open circles), (a) '^i^, 
(b) iijj, (c) (d) bji. Note that i?e(v['jj) = Im{Uj^^) = 7m($j^) = J?e(fejj) = 0. Because the Prandtl number is small, note 
rather sharp boundary layers appearing in ^l/^^ and iij^. The inset in (a) resolves the boundary layer behavior near x = —1. 



where T = e^t, Z = ez and the coefficients are defined in Appendix B. 

We stress here that the decoupling of these two equations is the result of translational (x-) symmetry of the thin- 
gap problem, but it cannot be guaranteed for a case in which, e.g., curvature terms have to be retained. Eq. (|19p is 
the diffusion equation and its physical implications are quite trivial. It indicates that the contribution of the above 
mentioned neutral mode to the azimuthal field perturbation will simply decay on a time-scale associated with the 
system's size and the smaller of either R„i or 1/CTZ - the meaning of the latter possibility will be explored in a 
forthcoming work. In contrast, equation is the well-studied real Ginzburg-Landau equation (see, e.g. [23l[2^[25j) 
which can exhibit non-trivial behavior in both the amplitude and phase of the envelope function A. The phase can 
lead to interesting dynamics emerging from Eckhaus-like instabilities, however in the present study we care only about 
the behavior of the amplitude's magnitude, i.e. the modulus of A. We shall thus agree henceforth to mean |^|, when 
writing A. Further discussion on phase dynamics can be found in the concluding section of this paper. 

A real amplitude A in the real GLE has two stable spatially uniform steady solutions, A{Z,T) = ztAg, and one 
possibly unstable solution, A = 0, as can be easily verified. Depending on the boundary conditions, the system 
typically relaxes to one of the steady solutions or, possibly, splits into two regions (the plus and minus values of Ag) 
with a front separating them (see e.g., [26,] for an example of a system of this kind). 

From (fT5|) it is apparent that the saturation amplitude is Ag — \f\fa. and thus its determination calls for the 
computation of the relevant coefficients. As discussed before this has to be done numerically. The details of this 
calculation are given in Appendix B and some representative results (for the parameter values g = 3/2 and C = 0.08) 
are displayed in Figure [3l Panel (a) demonstrates the weak dependence of Rm and Q, and of the coefficients A and 
-D, on Vra (for Vm ^ !)• In contrast, the coefficient a, whose numerical values are shown in panel (b), has a power 
law dependence on Vm in the same interval. Thus, the dependence of the saturation amplitude on Vm is essentially 
governed by a. The appropriate scaling for Vm ^ 1 is = A/a ~ 1/a (given the very weak dependence of A on 

Vm). 

The analysis sketched out in Appendix D shows that the dominant terms in the expression for a are such that 
a ^ VTa"^^ 1 for Vm ^ 1 This scaling fits very well the numerical results in Figure [5)3 (solid line). We thus obtain the 
following scaling behavior for the square of the saturation amplitude 

Al ~ V'^l^ (or ^2 - Vr^l'^ for fixed both for < 1. (20) 

The physical effects that these dominant terms are reflecting can be traced in the asymptotic analysis as resulting 
from the nonlinear radial advection of the second order azimuthal velocities Ux\dxUy2 and the creation of the azimuthal 
field due to the shearing of the radial perturbation field bxidxUy2- Note that in UMR06 we were able to obtain (from 
not fully consistent boundary conditions for this problem) the analytical result A"^ ~ Vm (or ^ Tl~^ for fixed i?„i). 
Thus we see that the implementation of more realistic boundary conditions that are appropriate for the thin-gap 
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FIG. 3: (Color online) Coefficients, parameters and saturation amplitude squared, as a function of Vm, for q — 3/2, C = 0.08. 
(a) Plots of Rm, Q, -D, A as a function of Vm- Note the weak sensitivity on Vm- (b) Values of a (diamonds) as a function of 
Vrn- Plotted as well is the functional dependence a ~ Vm*^^ (solid line) discussed in the text. 




FIG. 4: (Color online) Contours of the perturbation streamfunction (^P) and azimuthal velocity (uy) in the x — z plane, up 
to and including order e^. The vertical scale of the plots corresponds to one critical wavelength, 2tv/Q. The parameters are 
Vm ~ 0.001, C — 0.08, Rm ~ 4.9 with q = 3/2. We take e — 0.5 and the amplitude A = 0.07 (which is its saturation value for 
this case). Cuts along constant values of z are shown in the left panel (and dashed lines in the contour plots on the right ). 
Note that Uy is the velocity disturbance about the steady profile —qx. 
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Azimuthal B-field Contours b 
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Flux Function Contours <l> 




FIG. 5: (Color online) Same as Fig. |4]for the azimuthal field (top panels), by, and the flux function, <1> (bottom panels). 



niTC problem does not alter the general qualitative trend - saturation amplitude increasing with Vm (or decreasing 
with TZ for fixed Rm) - uncovered in UMR06, nor its implications. It merely alters (slightly) the power of this basic 
dependence. 

In Figure S] we plot the azimuthal velocity Uy{x, z) and the streamfunction "^{x, z) of the perturbation, calculated 
by our asymptotics to order . This has to be understood as the modification on top of the basic mTC configuration, 
which together constitute the steady saturated state. The presence of boundary layers near the channel walls is 
clearly apparent. In Appendix El we estimate that the boundary layer sizes scale as ~ Vrl^ and this is quantitatively 
consistent with the increase in power of the scaling from Aj. ^ Vm (as found in UMR06, where the boundary layers 
were essentially neglected) to Aj. ^ Vm^ here. The crucial ingredient in determining the scaling of Ag is, as we have 
seen, the scaling behavior of the coefficient a, which in turn is affected by the boundary layer width through its 
dependence on the relevant x-eigenf unctions (see Appendices A and D). 

In Figure [5] we display the perturbation's azimuthal field, hy, and its poloidal flux function, <i>(x, z), in a manner 
similar to the previous figure. Note that we do not see prominent boundary layers in the magnetic field perturbation; 
this is the result of the boundary conditions imposed ifT?)) . Whereas three velocity boundary conditions are imposed 
on each side (ensuring zero perturbation velocity on the boundary), only two such conditions on the magnetic field 
perturbation are enforced {hx = 0,dxby = 0). It is so because precisely ten conditions in all are required, otherwise 
the problem would be ill-posed. 

Finally, we turn to the evaluation of the angular momentum transport (a key question in assessing the MRI's role as 
the driver of accretion in astrophysical systems). The total (local) stress resulting from the perturbation is composed 
of the Reynolds and Maxwell stresses and, in our notation, has the form (see, e.g., Qjfsoj, UMR06) 

T{x, z) = Tr + Tm = UxUy - Cbxby. (21) 

As in UMR06, we may define a quantity measuring the average total angular momentum transport in the domain, 







/ dz 


J T{x,z)dx 


J-7T/Q 
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Total Stresses T(x,z) in domain 




-0.001 



FIG. 6: (Color online) Same as Fig. |4]for the total stress, T{x, z) (top panel), and the vertically integrated stress, ^ ~ J Ti^, z)dz 
(bottom panel). 



The quantity J can be thought of as analogous to similarly defined quantities used as a measure of the effective 
viscosity parameter due to dynamical fluctuations in active fluid media (for a recent purely hydrodynamic example, 
see [l3l) - be it either in a fully turbulent state or otherwise. In our problem, this quantity in the saturated state can 
be written to leading order as 



(22) 



where Jq = JoiRrm QjC, q) is of order unity. Given the behavior of the saturated amplitude in the V„i ^ 1 limit, it 
follows that the average angular momentum transport scales like 



(23) 



to leading order. Finally, we show in Figure [5] the distributed stress T{x,y) over the domain and the vertically 
integrated stress, defined by S(a:) = jTdz. 



V. DISCUSSION AND SUMMARY 



In this paper we have presented a full exposition of a weakly nonlinear asymptotic analysis of the MRI for a 
viscous and resistive flow in the thin-gap magnetic Taylor-Couette configuration. Our previous work (UMR06) 
employed mathematically expedient, but not fully consistent boundary conditions for this problem, so as to allow 
for transparent analytical evaluation of the envelope equation coefficients. Here we have used consistent and realistic 
boundary conditions for the mTC setup. As a result, the calculation is more involved. We have nevertheless found 
(as anticipated in UMR06) that in the thin gap limit the amplitude of the disturbances saturates at a value that 
decreases with decreasing magnetic Prandtl number, Vm- Moreover, the emergence of boundary layers actually makes 
the Pm-dependence of the saturation amplitude, and thus the average angular momentum transport more severe. 

Our results should be put in the proper context. They are valid close to instability threshold and in a confined 
system (mTC). Most previous studies of the MRI in the nonlinear regime (both numerical and analytical) followed 
the evolution of channel modes - exponentially growing, radially independent modes (see |5j]), which happen also to 
be exact solutions of the nonlinear equations for the perturbation, in the SB formulation under periodic boundary 



11 



conditions, i.e. in an open system. It is thus only natural that the channel modes have been identified as the dominant 
dynamics and their evolution perceived as a crucial ingredient in the nonlinear saturation of the instability. Goodman 
& Xu [2^ showed that the channel modes ultimately become unstable and break up. The asymptotic study of 
nonlinear saturation performed by Knobloch & Julien '271] was also based on a state dominated by channel modes. In 
these works, as well as the recent local modeling of MRI angular momentum transport [13, [Mj, results were compared 
with numerical simulations of an open SB (undoubtedly dominated by dynamics ari sing from the nonlinear evolution 
of channel modes). Note, however, that in the global approach of Kersale et al. 132', 'ss'l, the explicit inclusion of 
boundary conditions and curvature terms broke the radial symmetry of the problem (which is necessary for the channel 
modes to be manifested). These authors found numerically (using a spectral code) that the form of the saturated 
state critically depends on the boundary conditions adopted and, in any case, is not a "trivial" Keplerian state with 
developed MHD turbulence on top of it. 

From the vantage point of the linear theory followed here (as well as the SB investigations of the past), the MRI 
takes place primarily because the term supplying the tension, i.e. a perturbed azimuthal i3-field, arises from the 
sheared conversion (by the background flow) of a perturbed radial magnetic field, emanating from the bending of 
the background vertical field. The strength of the resulting destabilizing torque is related to the magnitude of q 
(measuring the local stretching) and the magnitude (squared) of the global vertical _B-field (representing that basic 
source of tension which is being stretched by the shear). Nonlinear saturation of a linear instability can generically 
be achieved by increased dissipation, by the modification of the linearly unstable base state so as to push it back to 
stability, or a combination of both. 

In the problem studied here, we have considered the marginal MRI mode (i.e. with growth rate 0), as a function of 
all free parameters, save q, which has been fixed to 3/2. We find that the saturated azimuthal velocity disturbance 
provides an effective positive radial gradient, q' > 0, through the bulk of the flow (see Figure S]). Thus the effective 
overall q in the saturated state is qcs = q — q' < 3/2. The magnitude of the effective gradient reflects the manner 
in which the modified gradient couples to the background field which is being stretched and is responsible for the 
instability. In our case, q' is positive and thus reduces the initial destabilizing shear, but not sufficiently to cancel it 
entirely. It has to be noted, however, that the saturated state is noi just the base flow with reduced shear. It includes 
also extra poloidal and azimuthal fleld, as well as poloidal velocity. This steady state is thus more complicated; the 
presence of velocity boundary layers complicates it even further. It is thus not trivial to identify a simple process for 
the saturation "mechanism" in this case. We note that our results share similarities with the saturation mechanism 
proposed by Knobloch & Julien [271 ^^r the saturated MRI state developed, in a particular asymptotic regime, from 
the unstable channel modes discussed above. 

We have followed into the weakly nonlinear regime a dissipative system, which was in a marginal balance and 
obtained a steady saturated state from a reduction of the shear, in places over the domain where it counts the most 
(in terms of azimuthal field production), and from the emergence of a steady flow and magnetic fleld conflguration. 
In terms of dissipation, it is instructive to consider the energy relationship In our steady saturated state the 

first integral is just oc / E{x)dx and therefore is positive (see the bottom of Fig. [6]). As dtE in this (steady) saturated 
state must be zero, the sum of the two dissipative integrals must be equal to the first one. We have verified that it is 
indeed so. 

We have not considered in this paper phase dynamics, which is an inherent feature of the more general envelope in 
complex GLE. Phase dynamics may be rich, in particular in two and three dimensions, admitting well-known pattern 
instabilities like Eckhaus and Zig-Zag and these, in turn, can lead to effects like phase turbulence and complicated 
defect dynamics 0, [2^. In what is considered here, where the coefficients of the one-dimensional GLE are real, all 
that remains of the above is just a possibility of an Eckhaus instability. This may merely introduce some non-steady 
readjustment to the overall pattern phase, but it leaves unaltered the overall amplitude scale of the basic pattern that 
emerges. In particular, our system is open in the z dimension and thus there should be no difficulty for the phase to 
adjust itself to a stable value (see [l^l p. 200). Because we are interested here in the scaling of the transport (which 
is expressed by an integral of the envelope over the domain), phase dynamics (although interesting in its own right) 
does not influence this measure and we have thus considered only the modulus of the envelope. 

Our results and flndings here should ultimately be compared to experiments and numerical simulations accompa- 
nying them. Extension of this type of analysis to a wide-gap mTC configuration is possible, but the results will be 
somewhat more complicated than those presented here, due to the inclusion of curvature terms. Preliminary calcu- 
lations indicate that the evolution of the perturbation amplitude in this case is governed by two coupled envelope 
equations (see Appendix B). The properties of the saturated state, however, appear similar in their salient features to 
the ones explored in this paper. The case of an initial helical field, for which experimental detection of the MRI has 
recently been reported [21| . can also be investigated in the weakly nonlinear asymptotic formalism employed here. It 
will the subject of future work. 

Further analytical investigations of the nonlinear MRI, of the kind reported here, will contribute toward assembling 
a deeper understanding of this important instability. Such investigations may also help in addressing the issues of the 
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effect of numerical resolution upon the resulting dynamics. In particular it could be useful to conduct simulations 
for, say, a fixed value of the magnetic Reynolds number (well below any contamination by numerical dissipation) and 
examine if and how does the transport change with resolution. Numerical studies of the MHD turbulent dynamo 
problem (e.g., [sS, 'ss']) have shown that such considerations are very important. The understanding of the role that 
the MRI plays in astrophysical disks, which in its full generality is a formidable problem, may be enriched by the 
experimental, analytical and numerical studies of simpler systems. 
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APPENDIX A: THE LINEAR SCALE OF THE BOUNDARY LAYER 

The best way to identify the scalings that are appropriate for the boundary layer is to rewrite (|16p as a single 
equation for, say, the streamfunction \E'. Setting the time-derivative to zero results in 

rvl/ = [{V,niDl - Q2)2 ^ cRlQ^f {Dl - Q^) + Rt2qCQ'' - ujIQ^DI - Q^fR^^] * = 0, (Al) 

where Wg = 2(2 — q) and where the simplifying notation = d/dx is also used. The operator is tenth order in 
derivatives. Inspection of its form suggests that retaining only the terms of (jAl[) that are dominant (for Vm ^ 1) in 
a small region of size with A > (the total x-domain size is 2 in our units) at either of the two boundaries gives 

{VlD'^ -ujIQ^DIrD^^Q. (A2) 

Treating all quantities as being of O (1) except for Vm, we can now see that the value of the exponent A must be 1/3. 
More explicitly, we consider a boundary layer by rescaling the x-coordinate around the boundaries at x — ±1. We 
define ^ = Vm^{x =F 1) and insert this into ()Aip revealing 

_ p-4A^2Q2^2^^4) ^ ^ ^ (p2-8A^ • • • , l) = 0, (A3) 

where = d/d£_. As Vm ^ 0, a distinguished balancing limit (see, e.g. [3^) may be achieved when 2 — lOA = — 4A, 
or when A = 1/3. In this case, all other terms in the boundary layer region are sub-dominant to the two terms 

1/3 

remaining. Thus, in the limit Vm <C 1, the size of the boundary layer scales as Vm ■ 



APPENDIX B: DERIVATION OF THE GINZBURG-LANDAU EQUATION 

To help facilitate the development of the weakly nonlinear theory we rewrite the equations of motion in the following 
way, 

VdtV -f N = £V -f e^gv, (Bl) 

in which 

£ = Lo + LiS, + Udl + Udl + Udt, (B2) 
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and where the matrices are defined as 

n 

Lo 



Li = 
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'^m '^x ' 
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dl 
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(B3) 



At marginahty, the background vertical field = \. The degree of linear instability is thus governed by the small 
parameter e defined by 

= 1 - Bo- (B4) 

The vectors above are defined by 

V= N= (iVvi„iV„,iV$,iVfc)T. 

We assume that during the nonlinear development there are two vertical scales emerging in the problem, namely z and 
Z = ez. We also assume that as one tunes the vertical field parameter into the MRI unstable state, the temporal 
response scales similarly to e^. Thus we say, for example for the streamfunction, that ^ = ^'(z, ez, e^t), and similarly 
for the other physical variables. If we assume that the solution forms follows V = V(a;, z, Z, T), then all operators in 
(|B1I) are re-expressed by applying the replacements 

d, — >d, + ed^ , dt — >e^d^. (B5) 

Thus, we now have 

t^Vd^Y + ^ = CV + eCxd^Y + e^C'^dl'V + e^QV + O (e^) , (B6) 

where 

Ci = U + 2L2d, + 3L3dl + AUdl (B7) 

£2 = 2L2 + 6L35^ + 12L49f. (B8) 

The boundary conditions, aside from periodicity in the vertical, are 

d^^ = 9^* = uy = d,^ = d^by = 0, at x = ±1. (B9) 

We expand all quantities in a perturbation series 

$ = e$i + e^$2 + e^'1'3 H 

* = e*i + e^vpa + 6^*3 ^ 

Uy — €Ul + e^U2 + e'ui, + ■ ■ ■ 
by = ebi + 6^62 + £^63 H 

or in other words 

V = eVi + e^Va + e^Vg + O (e^) . 
With the above expansion and multiple scaling Ansatz, it follows that the nonlinear terms are expressed in the series 



(BIO) 
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in which N = (iV'*^ N^'^\ iV^*', iV^''^)'^. In component by component form, these expressions are exphcitly given by 

(Bll) 
(B12) 
(B13) 
(B14) 









+ 








1- e^TVg^"' - 




[^') 


iV(*) 








{^') 


Nib) 











where 



TV^^*) = J(*i,V2*i) -CJ($i,V2$i) (B15) 



N^""' - J(*i,ui)-CJ($i,5i), (B16) 
4' 



Tvf^ = -J(<i>i,*i), (B17) 



^i'^ = J(«'i,6i)- J($i,ui), (B18) 



and 



TV;^*^ = 7(^-2, V^^i) + J(*i,V2*2) -CJ($2,V2$i) -CJ($i,V2$2) 

v2*i) -cj($i, v2$i) - 2j(*i,9,a^^'i) - 2j($i,a,a^$i), (bi9) 

TV;^"^ = J(*2,ui) + J(*i,u2)-CJ($2,6i)-CJ($i,62)-+J(*i,Mi)~CJ($i,6i), (B20) 
ivf) = -J($2,*i)- J($i,*2)- J($i,*i), (B21) 
TVf = 7(^-2, &i) + J(*i, 62) - J($2, ^^i) - «2) + &i) - J{^i,u,), (B22) 

in which J(/,g) = d^fdxQ - dxfd^g (remembering also that J{f,g) = dzfdxQ - dxfdzg). 
To O (e) we have 

£V, = 0. (B23) 

The solution to this is the marginal case investigated in the text. We write its general solution form as 

V, = A{T, Z)Vii e*'3^ + c.c. + B{T, Z)V,, , (B24) 

where 'V^^{x) = {'^^^,u^^,^^^,b^^)'^ , U^g = (0,0,0, 1)'^ and A and B (not to be confused with the magnetic field) are 
envelopes (amplitudes). Keeping in mind the above cited solution form, the boundary conditions at this order are 

= L'x*!! = = $11 = Dxb,, =0, at X = ±1, (B25) 

where = d/dx. Note that the constant azimuthal field symmetry discussed in the main text is embodied in the 
final term of l|B24[) . i.e. B{T, Z)l]-^-^. We also call attention to the fact that though all of these functions are order 1, 

they all (especially -u^J show the presence of boundary layers in a region O (Vm^^ close to the two boundaries. 

At order e^, the equations are 

£V, ==N2 -^i^^V,. (B26) 

The solution at this order is written as 

V, = A2V,,e'2Q. + a^AV,,e'«^ + lApV,, + d.BV,, + c.c. (B27) 

where = (^'2, "2, *2, ^2)"^ and in particular, Y^^ix) = (*22, U22, $22, ^22)'^, ^2ii^) = (*2i, "21, *2i, ^21)'^, 
Y^oix) = (^1/20, U20, ^"20, ^20) ,and 'iJ2o(^) = (0; '"20, 0, 0) . N2 contains no terms resonant with e*"^^ (see below) 
but because the expression ^CiC^^V^ does contain such a term, in order for there to be a bounded solution at this 
order with the required boundary condition, the following solvability condition (the vanishing of an inner product 
[33) must be satisfied: 



Vt • Cid^Y,) = §- I ! ''^ dxdzY^ ■ Cid^Y, = 0, 



(B28) 
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FIG. 7: (Color online) The adjoint solution Y"^ for C = 0.08, Vm = 0.001, q = 3/2 with ~ 4.9 and Q ~ 0.75. The open 
circles are the values on the Chebyshev grid while the solid lines are fits. 



where is the solution to the adjoint operation 

£^V^ ^ 0. (B29) 

The adjoint operator is given by 

£t = Lo - Ljd, + L2dj - L^'^d^, + Udj. (B30) 

and is so written since 

Lo = (Lo)'^, Li = (L2)'^, L4-(L3)'^. 

The adjoint solution = V^e**^^ + c.c., in which = (vj/t^ -^t ^t^T -g g^,,]-^ f^]^g^i satisfies the boundary 
conditions 

= 9^*^ = = $t = Q^i,i = 0, at X = ±1, (B31) 

in addition to periodicity in the vertical direction. In Figure we display an example of V^. We note that it is also 
an even function with respect to x. The solvability condition (|B28p is automatically satisfied on account of the choice 
of Q and Rm, as discussed in the main text. To complete this exposition, we explicitly write out the nonlinear terms 
appearing in N2, 

N2 = A^N22e'^^' + lApNao + c.c, 
where N22 ^ (7V(f ), 7V(f ), M*), iV(f N^o ^ (iV(*), iV(*), iV(*), M*))^ and 

-C [(*Q<i>,J • D^{Dl - Q^)^,, - (i?,$,J • iQ{Dl - g2)$^J 
Nlf = {iQ^,,)-D,b,,~{D,^,,)-iQb,,^C[{iQ<^,,)-D,u,,-{D.,<^,,)-iQu,,], (B32) 



and 



-C [(zQ$^J . [D^[Dl - g2)$,J* - (i^.<i>,J • {iQ{Dl - Q2)$^J*] 
^io^ = • {D^KT ~ {D.■^^^) ■ WK,r - • {D.U.J* ~ . itQu,,] 



(B33) 



16 




FIG. 8: (Color online) The second order solutions V22, V21, V20 for C = 0.08, Pm = 0.001, q = 3/2 with ~ 4.9 and Q ~ 0.75. 
Note the pronounced boundary layers in U22,tt2i, U20 and, to a lesser degree, ^'21. 



Superscript"*" on any given quantity (i.e. /*) denotes the complex conjugation of the said quantity. We solve for 
the quantities , V^j^ , V^j, using the Chebyshev collocation technique developed for the linear theory and show an 
example of these results in Figure [H We note that the functions of V„ , Y^g are odd with respect to x while those 

of V,, are even with respect to x. We also call attention to the fact that since iV^^^ + {Ni^)* = 0, is also zero. 

This is significant because it is really a consequence of a second condition that must be met to ensure the existence 
of a solution at this order. In particular, inspection of the z independent component of the equation describing the 
evolution of 62, i-e., 



(6) 



(B34) 



shows that in order for there to be a solution to b^^^ which satisfies the boundary conditions d^b^g = at x = ±1, a 
condition must be met with respect to the terms on the right-hand side of (jB34p . This criterion is most simply seen 
by (i) integrating this equation from x = — ltoa; = l, (ii) applying boundary conditions, (iii) leaving the requirement 



= 1^1^ 



N, 



(b) 



dx. 



However, this relationship is automatically satisfied on account of the fact that i?e(A^2Q-') = 0, as noted above. We 
complete the solution to this order by explicitly writing out the form for U20 which is associated with dzB - c.f. (|B27I) . 
The equation governing its structure is 



n~^dlu2o = -c. 



(B35) 



yielding the solution 



"20 



-C7^(2;2 - 1). 



(B36) 



Finally, at order t^, the equations are 

Vd^\, + N3 = CV, + £ia^V, + C^dlY^ + ^V^. (B37) 
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FIG. 9: (Color online) The non-linear function N31 for C = 0.08, Vm = 0.001, q = 3/2 with ~ 4.9 and Q ~ 0.75. 

The terms of the nonhnear functional N3 have the more detailed following expansion 

N3 = A^Naae*^^^^ + ANgae'^'^" + ^l^pNgie^^^ + Ad^BN^^^e'^' + ^*c)^AN3o + c.c. (B38) 

In order for there to be a solution at this order, the same solvability condition discussed earlier must also be satisfied 
for this equation. Those terms in the above expression subject to the solvability condition are the ones resonant with 
gjQz ^ggg i^j^g definition of above). Our goal in this work is to satisfy this solvability condition, at this order. 
In turn, this means that the only two terms from the nonlinear expression that will concern us here with this first 
solvability (see below) will be the expressions involving N31 and N31. The explicit forms for these expressions are 
given in Appendix C (in order not to clutter this exposition). The main feature to note about these functions is that 
N31 is even with respect to x while N31 is odd. The ramifications of these facts are explained below. The solvability 
condition means taking the inner product of (jB37p with V^, revealing 

ad^A + cA\A\^ +dAd^B ^bA + hdlA, (B39) 

where 

a = <Vt . VY*J , c = (Vt . N5i) , S = (Vt . N*i) , b={Y^- (gV,J*) , 

h ^ (vt.f£iV,,+£2V,,)*\. (B40) 



Inspecting the expressions comprising N31 reveals that they are odd symmetric with respect to x. This result means 
that the expression c — and, consequently, it means that the evolution equation for the field A evolves according to 

d^A^ XA + DdlA- aA\A\'^, where X = b/a, D = h/a, a = c/a, (B41) 

independent of the second field quantity B. This decoupling is a direct consequence of the symmetries preserved in 
the thin gap Taylor- Couette limit. This will not be the case in the general-gap magnetized Taylor- Couette problem. 

As noted earlier, there is a second solvability condition, which must be enforced on the order e'^ _B-equation. This 
becomes necessary on the component of the solution for which there is no explicit z dependence. More explicitly, 
writing out that component of the equation: 



R-^d%o = d^B - R-^dlB~ u2odlB 



A*d,AN^, 



(b) 



(B42) 



The solvability condition can be readily inferred by integrating (IB42P from a; = — ltoa; = l and requiring that 
dxb-io = at x = ±1. This procedure is equivalent to taking the inner product of (|B37p multiplied by U^^, thus 
revealing the second envelope equation 

d,B=^-^+'^-^yiB~pdM\ 

where p = (^Re (j^-jo'^ ^ . However, the odd symmetry property of A^3q'' means that p = 0. Thus, the second solvability 
condition yields the simple diffusion equation, 

9TB=(^ + ^-^)dlB. (B43) 
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APPENDIX C: THE TERMS FOR N31 

We define for notatfonal convenience the symbols Vf = D^ — and V2 = -D^ — AQ"^. It follows that the expressions 



are 



and 



-C {(iQ$,J* • (i?.V2cI>,J + (zQa>,J . - (i?.<i>iJ* • {t2Q\/l^,,)} , (CI) 

iV(^-) = {t2Q^.,.,){D,u,,r -[{D^^,,)-{tQu,,r + {D.,'i',,)-{zQu,,)] 

-C{{tQ^,,y ■ (D^b,,) + (zQ*,,) • D^b,„ - (D^^,,)* ■ {i2Qb,,)} , (C2) 
ATf ) = - {iz2Q^,,)iD,^,,r - m^,,) ■ (*Q*^J* + . (iQu^,,)]} 

{(*Q*,J* • p.*,,) + {iQ^.J ■ D,^,„ - {D,^,,r ■ (*2Q*^J} , (C3) 
TVW = iz2Q^,JiD,b,,r-[iD,^,J-izQb,,r + iD,^J-izQb,J] 

- {{i2Q^,,){D^u,,r - • (iQuiJ* + {D,^,,) ■ Wu,,)]} 

- {W$,J* • {D^u,,) + (tQ^^JJJ^- • (*2Qu,j} (C4) 

iV^*) = (C5) 

iVir^ = W^I,^ J • i^xit.o (C6) 

TVff) = (C7) 

ivW = ^{tQ^,J-D,u,„ (C8) 



APPENDIX D: ON THE O ( P^''^^ ) DEPENDENCE OF a. 



There are a number of terms comprising the integral expression leading to the quantity a. Of these, there are a 
few that dominate its expression when Vm is small. In the following, we will sketch out one way to understand the 

O (Vm'^^^^ scaling behavior of a. (Note that because we consider the behavior of this system by holding TZm = Km 

fixed, we will speak about the general scaling dependencies of quantities on Vm and TZ~^ interchangeably as they are, 
in effect, equivalent under this constraint.) 

We noted earlier that the lowest order functions comprising V^j exhibit boundary layers of spatial extent O (Vm^^ 
for Vm ^ 1. Especially acute in this respect is the function for the azimuthal velocity perturbation, u^^. At the next 
order, we find that the equation for u^^ (i.e. the second order azimuthal velocity function with no dependence on the 
z coordinate) is simply 

n"d^,u,,=N^:l (Dl) 

Inspection of (|B33p shows that N^^'> is dominated by the underlined term containing the expression D^u^^. (Note 
that it is true that the other quantities also, in principle, have boundary layers as well but the one associated with 
DxU-^^ dominates - inspection of Fig. [5] readily shows.) This means that in the boundary layer regions, N^'^^ scales as 

O {Vm"^^^ while in the interior it remains O (1). It follows from inspection of (|Dip that u^^ scales as O {Vm^^^^ in 

the boundary layers and O {Vm^) in the bulk interior. 

Now wc turn to an inspection of the expression leading to a, namely term c of (jB40p which is composed, in part, of 
the integrals over the domain of the products ■ N^^^ and b^ ■ N^^'^K Since and b^ remain O (1) over the entirety 
of the domain, it remains for us to evaluate the behavior of iV^^"-* and N^'^'' over the domain. Although there are 
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several terms that contribute, one of the terms is most dominant: the underhned expressions in (jC2[) . which involves 

— 1/3 

DxU^g. Because there is a derivative, the scale of D^u^g gets amplified by another factor of Vm in the boundary 
layer region. Given what we have established thus far about the character and profile of w^^, it follows that D^u^g is 

O (Vrn'^'^^ in the boundary layer regions while it is O (Pm^) in the interior. It means, therefore, that the profiles of 

iVjj^"-' and N^^'^^ similarly reflect this character on the x domain. 

Thus, for instance, to determine the order of magnitude of the integral of b^N^'^'> over the domain, we should break 

1/3 

up the integral into parts separating out the interior region and boundary layers. Writing S = Vm we have 

b^N^^'^dx= b^N^^Ux+ b^Ni^^Ux+ I b^N^^Ux 



1 J-l J-l+S Jl-5 



In other words, because the length scale of the interior zone is O (1) and the scale of b^N^^^'^ = O (V,-^) in that 
region, the contribution to the integral from this part is O ("Pm^)- On the other hand, because the length scale of 
the boundary layer(s) is O (Vm^^ , while the value of ^t^y^W scales as O (Vm^^^^ in those regions, it follows that 

the contribution to the total integral from these zones is O (Vm^^^^ ■ The same reasoning follows for the integral of 
over the domain. 

The nonlinear readjustment occurring in the boundary layers dominates the scale of a and we can conclude that the 
dominant process leading to saturation occurs in the boundary layers. We note also that had there been no boundary 
layers, then a would scale as O {'P~^) because of the scale of u^^, which is always at least O iV^) on account of (|D1|) . 
This directly relates to the problem investigated in UMR06, in which boundary layers are suppressed on account of 
the boundary conditions employed in that study. In that case, the saturation process gets contributions from the 
entirety of the domain and not just the boundary layers. 
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